#!/usr/bin/perl


use strict;
use Getopt::Long;
use FindBin qw($Bin);
use lib "$Bin";
use lib "$Bin/../lib";
use File::Basename;
die "usage:perl $0 <phy> \n" unless @ARGV == 1;

my    $data       = "dna";
my    $rates      = "gamma";
my    $nst        = 6;
my    $ngen       = 100000;
my    $samplefreq = 100;
my    $aamodel    = "Jones";
my $burnin = ($ngen / $samplefreq) * 0.25;
open(IN,  "$ARGV[0]")          or die "Can't open :$ARGV[0]\n";
open(OUT, ">input.nex") or die "Can't open :input.nex\n";
my $firstline = <IN>;
my ($species_number, $sequence_length);
if ($firstline =~ /^\s*(\d+)\s+(\d+)/) {
	$species_number = $1, $sequence_length = $2;
	print OUT
	"#NEXUS\n\nbegin data;\n\tdimensions ntax =$species_number nchar=$sequence_length;";
}
print OUT "\n\tformat datatype=$data interleave =no missing=? gap=-;\n\tmatrix\n";
while (my $eachline = <IN>) {
	print OUT "\t$eachline";
}
print OUT "\t;\nend;\n";
print OUT "begin mrbayes;\n";
print OUT "\tset autoclose=yes nowarn=yes;\n";
if (lc($data) eq "dna") {
	print OUT "\tlset nst=$nst rates=$rates;\n";
}
elsif (lc($data) eq "protein") {
	print OUT "\tprset aamodelpr=fixed($aamodel);\n\tlset rates=$rates;\n";
}
print OUT
"\tmcmc ngen=$ngen printfreq=100 samplefreq=$samplefreq nchains=4 savebrlens=yes;\n";
print OUT
"\tsumt filename=input.nex contype=allcompat burnin=$burnin;\n";
print OUT "\tquit;\nend;\n";
close OUT;
print STDERR "/mnt/sdb/zengl/bin/rna-pipe/Compare-rna/ploy/mb -i input.nex\n";
